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INTRODUCTION 


The  direct  application  of  the  finite  element  method  to  crack 
problems  was  studied  by  a number  of  investigators  [1-3].  No  special 
attention  was  given  to  the  singular  nature  of  stress  and  strain  of  a 
crack  tip.  Because  of  the  large  strain  gradients  in  the  vicinity  of  a 
crack  tip,  it  requires  the  use  of  an  extremely  fine  element  grid  near 
the  crack  tip.  By  comparing  the  finite  element  result  of  displacement 
components  or  stress  components  at  a nodal  point  with  the  corresponding 
asymptotic  result  of  displacement  or  stress  components  at  that  node,  the 
stress  intensity  factor  can  be  estimated.  The  estimated  values  of  a 
stress  intensity  factor  vary  over  a considerable  range,  depending  on 
which  node  is  taken  for  computation.  This  results  in  poor  estimates  if 
displacements  are  taken  at  nodal  points  either  very  close  to  or  far 
away  from  the  crack  tip. 

An  improved  finite  element  technique  was  developed  by  Wilson  [4]. 

It  combined  the  asymptotic  expansion  of  displacements  in  a small 
circular  core  region  surrounding  a crack  tip  and  the  finite  element 
approximation  outside  a polygon  approximating  the  circular  arc  of  the 

^Swedlow,  J.  L.,  Williams,  M.  L.,  and  Yang,  W.  H.,  "Elasto-Plastic 
Stresses  and  Strains  in  Cracked  Plates,"  Proceedings  First  International 
Conference  on  Fracture,  1,  p.  259,  1966. 

^Kobayashi,  A.  S.,  Maiden,  D.  E.  and  Simon,  B.  J.,  "Application  of  the 
Method  of  Finite  Element  Analysis  to  Two-Dimensional  Problems  in 
Fracture  Mechanics,"  ASME  69-WA/PVP-12  (1969). 

3Chan,  S.  K.,  Tuba,  I.  S.  and  Wilson,  W.  K. , "On  Finite  Element  Method 
in  Linear  Fracture  Mechanics,"  Engineering  Fracture  Mechanics,  2,  p.  1, 
1970. 

^Wilson,  W.  K. , "Combined  Mode  Fracture  Mechanics,"  Ph.D.  Dissertation, 
University  of  Pittsburgh,  1969. 
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core  region.  The  displacement  fields  obtained  from  these  two  approx- 
imations are  not,  in  general,  continuous  along  the  asymptotic  expansion- 
finite  element  interface  except  at  discrete  nodal  points. 

An  alternative  finite  element  approach  to  crack  problems  is  the 
use  of  special  elements  in  the  region  of  the  crack  tip,  e.g.  [5-7]. 

In  [5],  Tracey  employs  quadrilateral  isoparametric  elements  which 
become  triangular  around  the  crack  tip.  The  displacement  functions  of 
the  two  types  of  elements  are  selected  such  that  displacements  are 
continuous  everywhere,  and  the  near  tip  displacements  are  proportional 
to  the  square  root  of  the  distance  from  the  crack  tip. 

Henshell  and  Shaw  [8]  and  Barsoum  [9]  showed  that  special  crack 
tip  elements  were  unnecessary.  For  two-dimensional  8-node  quadrilateral 
elements,  the  inverse  square  root  singularity  of  the  strain  field  at 
the  crack  tip  is  obtained  by  collapsing  quadrilateral  elements  into 
triangular  elements  and  placing  the  mid-side  nodes  at  quarter  points 


5Tracey,  D.  M. , "Finite  Elements  for  Determination  of  Crack  Tip  Elastic 
Stress  Intensity  Factors,"  Engineering  Fracture  Mechanics,  Vol . 3,  1971. 

6Blackburn,  W.  S.,  "Calculation  of  Stress  Intensity  Factors  at  Crack 
Tips  Using  Special  Finite  Elements,"  The  Mathematics  of  Finite  Elements 
and  Applications,  Brunei  University,  1973. 

7Benzley,  S.  E.  and  Beisinger,  A.  E.,  "Chiles-  A Finite  Element  Computer 
Program  that  Calculates  the  Intensities  of  Linear  Elastic  Singularities," 
Sandia  Laboratories,  Technical  Report  SLA-73-0894,  1973. 

8Henshell,  R.  D.,  and  Shaw,  K.  G.,  "Crack  Tip  Finite  Elements  Are 
Unnecessary,"  International  Journal  for  Numerica^  Methods  in  Engineering, 
Vol.  10,  1975. 

9Barsoum,  R.  S.,  "On  the  Use  of  Isoparametric  Finite  Elements  in  Linear 
Fracture  Mechanics,"  International  Journal  for  Numerical  Methods  in 
Engineering,  Vol.  10,  1976. 
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from  the  tip.  The  quarter-point  quadratic  isoparametric  elements,  as 
singular  elements  for  crack  problems,  have  been  implemented  in  NASTRAN 
by  Hussain  et  al  [10]. 

In  order  to  reduce  the  computer  core  requirement  and  to  simplify 
the  modeling  of  a structure,  better  known  but  lower  order  finite 
elements  have  been  abandoned  in  favor  of  cubic  12-node,  isoparametric, 
quadrilateral  elements  as  described  by  Zienkiewicz  [11].  In  this  paper, 
the  concept  of  quarter-point,  quadratic,  isoparametric  elements  is 
extended  to  12-node  cubic  isoparametric  elements.  The  correct  order  of 
strain  singularity  at  the  crack  tip  is  achieved  in  a simple  manner  by 
collapsing  the  quadrilateral  elements  into  triangular  elements  and  by 
placing  the  two  middle  nodes  of  a side  at  1/9  and  4/9  of  the  length  of 
the  side  from  the  tip.  The  12-node,  isoparametric  elements  have  been 
implemented  in  NASTRAN.  Both  mode  I and  mixed  mode  crack  problems  are 
computed  by  NASTRAN  using  the  collapsed  elements  to  assess  the  accuracy. 
The  stability  of  results  is  discussed  when  the  collapsed  triangular 
elements  are  used. 

THE  1 2-NODE  QUADRILATERAL  ISOPARAMETRIC  ELEMENT 

A typical  12-node,  quadrilateral  element  in  Cartesian  coordinates 
(x,y)  which  is  mapped  to  a square  in  the  curvilinear  space  (£,ri)  with 
vertices  at  (±  1,  ± 1)  is  shown  in  Figure  1.  The  assumption  for 
displacement  components  takes  the  form: 

^Hussain,  M.  A.,  Lorensen,  W.  E.,  and  Pflegl,  G. , "The  Quarter-Point 
Quadratic  Isoparametric  Element  As  a Singular  Element  for  Crack 
Problems,"  NASA  TM-X-3428,  1976,  p.  419. 

^Zienkiewicz,  0.  0.,  The  Finite  Element  Method  in  Engineering  Science, 
McGraw  Hill,  London,  1971. 
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12 

u = I Ms,n)u- 

1=1  v 

/ 1 

12 

v=  l Ni(5,n)v1 

i=l  ' 

where  u,v  are  x,y  components  of  displacement  of  a point  whose  natural 
coordinates  are  £»n»  are  displacement  components  of  node  i and 

N^(?,n)  is  the  shape  function  which  is  given  by  [11] 

Ni(?’n)  = ^6  0 + 5^)0  + nn-i )C-10  + 9(e  + n2)][-io  + 9(s?  + nf )3 

+ (1  + )(1  + - n2 ) ( 1 - n2) 

+ (1  + Titled  + 9^i)(l  - ?2)d  - 5?)  (2) 


for  node  i whose  Cartesian  and  curvilinear  coordinates  are  (x.j,y.j)  and 
(5-jjn-j)  respectively.  The  details  of  the  shape  functions  and  the 
numbering  sequence  are  given  in  Figure  1. 

The  same  shape  functions  are  used  for  the  transformation  of 
coordinates,  hence  the  name  isoparametric, 

12  \ 

X = l M§,ti)x. 
i=l  1 

> (3) 

12  ( 
y = I N-(5,n)y. 

i=l  1 1 ; 


^Zienkiewicz,  0.0.,  The  Finite  Element  Method  in  Engineering  Science, 
McGraw  Hill,  London,  1971. 
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The  element  stiffness  matrix  is  found  in  the  usual  way  and  is 


given  by  [9,10] 


1 1 

M-J  / 

-l  -l 


[B]T[D][B]  det  |J|d£dri 


(4) 


where  [B]  is  a matrix  relating  joint  displacements  to  strain  field 


[B]  - [ • • • B.j 


9Ni 

9x 

0 

9N-: 


9Ni 

9y 

9N4 


9y  9x 


(5a) 


and  [D]  is  the  material  stiffness  matrix  and  is  given  for  the  case  of 


plane  stress  by 


1 v 0 

v 1 0 

0 0 (1  - v)/2 


(5b) 


in  which  E is  Young's  modulus  and  v is  Poisson's  ratio. 
The  Jacobian  matrix  [J]  is  given  by 


^Barsoum,  R.  S.,  "On  the  Use  of  Isoparametric  Finite  Elements  in  Linear 
Fracture  Mechanics,"  International  Journal  for  Numerical  Methods  in 
Engineering,  Vol . 10,  1976. 

^Hussain,  M.  A.,  Lorensen,  W.  E.,  and  Pflegl,  G.,  "The  Quarter-Point 
Quadratic  Isoparametric  Element  as  a Singular  Element  for  Crack 
Problems,"  NASA  TM-X-3428,  1976,  p.  419. 
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1 — 1 

C-. 

1 1 

II 

9x  9^ 

3?  3? 

9x  9y 

9ri  9n 

It 

...  9N-j  ... 

9T 

• ••  9N  ^ ••• 

3n 

1 

• • X • 

_ll 

• • • • 

—1. 

1 

L_  J 

whenever  the  determinant  of  [J]  is  zero,  the  stresses  and  strains  become 
singular  [8-10].  The  derivatives  of  shape  functions  are 

— + nn-j)[-10  + 9(5?  + nf)](-105-j  + 18?  + 27^5*  + 

9£  CO 0 1 III 

+ 5,(1  + 9nn-j ) (1  - n2)(l  - nf) 

* (i  * nrl1)(i  - 5*1(95,  - 2?  - 275152)  (7a) 

3N  * i 

2^  (1  + + 9<52  * n2)](-10n,  + 18n  + 27n^n2  + 9ni52 ) 

+ ^n,(l  + 9«i)0  - 52)0  - 5f) 

+ ^ 0 + 55,1(1  - n|)(9n,  - 2n  - 27n,n2)  (7b) 

THE  CRACK  TIP  ELEMENT 

In  an  8-node  quadratic  isoparametric  element,  Henshell  and  Shaw  [8] 
and  Barsoum  [9]  found  independently  that  the  strain  became  singular  at 
the  corner  node  if  the  mid-side  nodes  were  placed  at  the  quarter  points 


8Henshel 1 , R.  D.,  and  Shaw,  K.  G.,  "Crack  Tip  Finite  Elements  Are 
Unnecessary,"  International  Journal  for  Numerical  Methods  in 
Engineering,  Vol . 9,  1975. 

^Barsoum,  R.  S.,  "On  the  Use  of  Isoparametric  Finite  Elements  in  Linear 
Fracture  Mechanics,"  International  Journal  for  Numerical  Methods  in 
Engineering,  Vol.  10,  1976. 

^Hussain,  M.  A.,  Lorensen,  W.  E. , and  Pflegl,  G.,  "The  Quarter-Point 
Quadratic  Isoparametric  Element  as  a Singular  Element  for  Crack 
Problems,"  NASA  TM-X-3428,  1976,  p.  419. 
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of  the  sides  from  the  corner  node.  This  singularity  is  achieved  in  a 
similar  way  for  a 12-node  isoparametric  element  by  placing  the  two 
middle  nodes  at  the  1/9  and  4/9  of  the  length  of  the  sides  from  the 
common  node  of  two  sides. 

For  simplicity,  let  us  consider  the  singularity  along  the  side 
n = -1  of  Figure  1.  In  general,  the  cubic  mapping  functions  are 

(8) 

(9) 


x = aQ  + a^  + a2S2  + a3r 


u = b0  + M + b2^  + b3?2 


For  £ = -1,  -1/3,  1/3  and  1,  the  corresponding  values  of  x and  u 
are 

x = 0,  a£,  6£»  £ 

U = U-j  , u2,  u3,  u4 

The  constants  a's  and  b's  in  terms  of  these  values  of  x and  u are 


X,  , - 

a0  " 16  (-1 

+ 

9a  + 96) 

» 91 

% 

= 16 

(-1  - 27a  + 276) 

1 

) do) 

9£ 

a2  " 16 

(1 

- a - 6) 

, a3 

9£ 

" 16 

(1  + 3a  - 36) 

J 

_ 1 
16 

(-u1  + 9u2 

+ 

9u3  - u4) 

’ bl 

_ 1 
16 

(u^  - 27u2  + 27u3  - 

u4} 

(ID 

9 

16 

(u]  - u2  - 

U3 

+ U4) 

* b3 

9 

16 

(-u.|  + 3u2  - 3u3  + 

u4} 

To  have  singular  strain  at  x = 0 (£  = -1),  the  reduced  Jacobian,  ^ , 
must  vanish  at  £ = »1.  From  (8)  we  have 

= a,  + 2a9£  + 3a-S2  (12) 

d£  1 2 3 
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(13) 


For  E,  = -1 , = 0 leads  to  the  equation 

3 = 2a  + | 

. n . n dll 

In  order  to  have  the  inverse  square  root  singularity  for  ^ , 

^=if=<bi+2^+3b^>/i 


x must  be  a quadratic  function  of  E,  so  that  the  inverse  gives  ^ as  a 
function  of  x?2.  This  leads  to  a3  = 0 or 


1 + 3a  - 38  = 0 

(14) 

The  solution  of  (13)  and  (14)  gives 

a = 1/9  and  8 = 4/9 

(15) 

Equations  (8)  and 

(9)  become 

x = | (1  + 0*  or  5 « -1  + 2/1 

(16) 

u = ui  + \ (-1 1 Ui 

+ 18u2  - 9u3  + 2u4)Jf  + f (2u1  - 5u2  + 4u3  - 

u4}  I 

+ j (-u-|  + 3u2  - 3u3  + u4)(j)3/2 

(17) 

From  (17)  it  is  clear  ^ has  singularity  of  the  order  1 at  x = 0. 

The  inverse  square  root  singularity  at  x = 0 along  any  other  ray 
emanating  from  node  1 can  be  achieved  by  degenerating  the  quadrilateral 
element  into  a triangular  element  with  the  side  10,  11,  12,  1 collapsed 
to  a point  at  the  crack  tip  and  placing  grid  points  2,  9 at  1/9  and  3, 

8 at  M/9  from  the  tip  (Figure  2),  where  % is  the  length  of  the  sides 
corresponding  to  ri  = ± 1 . The  Cartesian  coordinates  of  nodal  points  are 
shown  in  Figure  2 . Using  (3), 
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x = f (1  + 5)2 f(n>a,3) 
y = f (1  + 5)zf 1 (n>a,3) 

where  f and  f'  are  abbreviations 

f(n.a,B)  = (1  - n)cosg  + (1  + n)cosa 
f'(n,a,8)  = (1  - n)sing  + (1  + n)sina 
The  Jacobian  matrix  is  given  by 


3x  3y_ 
35  35 

| (1  + 5)f(n,a,0)  | (1  + 5)f'(n,a,3) 

1 

9?|S? 

l 

4 (1  + 5)2 (cosa  - cos0)  4 (1  + 5)2(sina  - sing) 

8 o _ 

and  the  determinant 

1 0 1 = (1  + ?)3sin(a  - 8)  (19) 


This  shows  the  strain  is  singular  at  x = 0 (£  ■ -1)  along  any  ray  from 
x = 0 since  |J|  = 0 at  £ = -1  for  all  n« 

For  the  inverse  functions,  we  have 


li  3n 
3x  3x 

1L  In 
_ 3y  3y  _ 

2(sina  - sing) 

-2(cosa  - cosB) 

[j]'1  = 

= 

5,(1  + 5)sin(a  - 8) 

-4f ' (n,a,B) 

5,(1  + 5)sin(a  - B) 

4f(n>a>B) 

_5,(1  + 5)2sin(a  - B)  A(1  + 5)2sin(a  - B)  J 

(20) 


In  terms  of  polar  coordinates  (r,0)  we  have  from  (18) 

1+5=2^  , R = (rM)cos(e  - ^)/cos(^) 
n = tan(0  - 2y^-)/tan(^y^-) 


(21) 
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The  displacements  components  u,v  at  a point  (£,n)  of  the  trian- 
gular element  of  Figure  2 are 
12 

o - J,  N-j(?>ri)u-j  - Aq ( ri , u^ ) + Ai (n>Ui*)(1  + £)  ^ 


+ A2(n,ui)(l  + ?)2  + A3(n,u1.)(l  + ?)3 
12 

v = J Ni(?,n)vi  = Ag(n»v^ ) + A-, (n.v^ ) (1  + 5) 
+ A2(n5vi)(l  + £)2  + A3(n,v.j)(l  + 5)3 


(22) 


where 


AgCn.u^  = {2 (— 1 + 9n2)[(l  - n)u1  + (1  + n)u10] 

+ 18(1  - n2 ) [ ( 1 + 3n)un  + (1  - 3n)u12]}/32  (23) 

The  displacement  derivatives  are 

3A 

3u  = 3u3i+3u3ri=  1 -4f‘  (n,ot,B)  0 

3x  3?  3x  3n  3x  (1  + ?)2  £sin(a  - 6)  3n 

1 2 3Ai 

’ TTTi)  iVi'n(a  - 6){(Slna  ' Sl"B)fl'  ' 2f'(r'-°'6)  8T>  + "•  <24> 

3u=3u3i+3u3n=  1 4f(n,a»3)  ^0 

3y  3?  3y  3n  3y  (1  + ?)2  £sin(a  - 6)  3n 

3A 

77  j ^ T~~r -r  {-(cosa  - cos3)A,  + 2f(n,a, 3)~->  + •••  (25) 

(1  + £)  £sin(a  - 8)  1 3n 


where 

— — — = {(2  + 36n  - 54n2)u-i  - (2  - 36n  - 54n2)uln 
3ri  1 IU 

+ 18(3  - 2n  - 9n2 )u-| -j  - 18(3  + 2n  - 9n2)u]2}/32  (26) 

Similar  expressions  for  3v/3x  and  3v/3y  can  be  obtained  by  replacing 
u-j  with  v^. 
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It  can  be  seen  that  the  strain  field  is  singular  at  r = 0.  The 
order  of  singularity  is  (1/r).  The  leading  term  vanishes  together  with 


SAg/Sri  for  all  n if 


U1  ' U10  U11  U12 


V = V = V = V 

1 10  11  12 


(27) 


In  this  case,  the  order  of  singularity  becomes  (l//r).  This  is  analo- 
gous to  the  constraints  given  in  [12]  for  quadratic,  isoparametric 
elements.  Hence  we  have  two  types  of  strain  singularity  at  our 
disposition.  If  nodes  1,  10,  11,  12,  which  are  collapsed  into  one 
point  at  the  crack  tip  (Fig.  2),  are  tied  together  during  deformation, 
the  elastic  singularity  is  obtained.  If  these  nodes  are  allowed  to 
slide  with  respect  to  each  other,  then  the  strain  singularity  is  of  the 
order  of  (1/r),  the  perfect  plastic  singularity. 

Using  the  multiple  constraint  conditions,  equations  (27),  the 
displacement  components  at  (£,n)  relative  to  the  tip  may  be  written  in 
the  form 

U = lV  v/R[36Fi(n,u.j)  + F^ruu-j)  + 36{F2(n,u.j ) - F-|  (n,u-j ) }/r 

- 36F2(n,ui)R]  (28) 

v = lV  v'7r[36F-|  (n»v.j ) + F3(n»v1-)  + 36{F2(ri,v.j ) - F-|(n»v^)}/R 

- 36F2(n,vi)R]  (29) 
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where 


F-j (n»ui ) = (1  - n)(2u2  - u3)  + (1  + n)(2u9  - ug) 

F2(ri»Li-j)  = (1  - n)(-3u2  + 3u3  - u^)  + (1  + n)(-3ug  + 3ug 

F3(n,ui)  = 9(1  - n2)[(l  - 3n)u5  + (1  + 3n)ug] 

- (1  - 9n2)[(l  - n)u4  + (l  + n)u7] 
and  Fi(n»v*)  etc.  are  obtained  by  replacing  with  v^. 


DETERMINATION  OF  STRESS  INTENSITY  FACTORS 

The  collapsed,  triangular  elements  around  the  crack  tip  have  the 
correct  elastic  singularity  at  the  tip  if  all  nodes  at  the  tip  are  tied 
together  during  deformation.  Only  one  set  of  displacement  functions  is 
used  for  regular  quadrilateral  elements  and  the  collapsed  triangular 
elements.  Hence  the  continuity  of  displacement  components  is  insured 
throughout  the  region.  The  nodal  displacements  obtained  from  the  finite 
element  method  using  higher  order  polynomials  for  the  displacement  field 
should  be  quite  accurate.  In  this  section,  we  will  briefly  discuss 
various  techniques  to  estimate  the  stress  intensity  factors  from  the 
nodal  displacements  thus  obtained.  The  discussion  is  limited  to  the 
use  of  nodal  displacements.  Other  techniques,  such  as  J-integral,  the 
strain  energy  release  rate  etc.  are  not  included. 

The  well  known  classical  near  crack  tip  displacements  are  given 
by  [13] 


^Williams,  M.  L.,  "On  the  Stress  Distribution  at  the  Base  of  a 
Stationary  Crack,"  Journal  of  Applied  Mechanics,  Vol . 24,  1957. 
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n-1  n-1/2 

2Gu (e ) = l {(-1)  r [d  D n(n,e)  + a A (n,e)] 

n=l ,2, . . . 2n-l  ul  2n-l  ul 

+ (“1)  r td2nDu2^n’0^  + a2nAu2^n,0^}  (31) 

n-1  n-1/2 

2Gv(e ) = £ {(-1)  r [d  D (n,0)  + a A (n,e)] 

n=l,2,...  2n-l  vl  2n-l  vl 


n n 

+ (-1)  r [d  D (n,e)  + a A (n,0)]}  (32) 

2n  v2  2n  v2 

where 

Dui(n,0)  = (n  - l/2)cos(n  - |)e  - (k  + n - |)cos(n  - l)e  N 


*^u2 (n ’ 6 ) = ncos(n  " 2)0  - (k  + n + 1 )cosn0 
Aul(n,0)  = (n  - l/2)sin(n  - |)0  - (k  + n + l/2)sin(n  - 1/2)0 


Au2(n»0)  = nsin(n  - 2)0  - (<  + n - l)sin0  (33) 

D^^(n,0)  = -(n  - l/2)sin(n  - |-)0  - (<  - n + |-)sin(n  - ^-)0 


(n > 0 ) = -nsin(n  - 2)0  - (k  - n - l)sin  n0 
1 5 

Avl(n,0)  = (n  - 2-)cos(n  - jjOe  + (k  - n - l/2)cos(n  - 1/2)0 


Av2(n,0)  = ncos(n  - 2)0  + (k  - n + l)cos  n0  / 

in  which 

{(3  - v)/(l  + v)  for  plane  stress 

(34) 

3 - 4v  for  plane  strain 
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The  coefficients  d's  and  a's  are  to  be  determined  from  the  boundary 
conditions  of  a problem.  The  stress  intensity  factors  K-j  and  K2  are 
related  to  d-j  and  a-j  by 

K]  = -d-jv^F  , K2  = -a-j^TT  (35) 

In  (31)  and  (32), u and  v are  displacement  components  referring  to 
the  local  Cartesian  coordinates  with  crack  tip  as  the  origin  and  the 
crack  on  the  negative  x-axis.  A simple  transformation  can  be  used  to 
change  the  nodal  displacements  in  global  coordinates,  obtained  from 
the  finite  element  method,  to  displacements  in  local  coordinates.  For 
simplicity,  we  will  assume  the  local  coordinates  and  the  global  coor- 
dinates are  the  same. 

1 . One  Term  Expansion: 

If  we  retain  only  the  /r  term  and  drop  all  higher  order  terms 

in  the  right  hand  sides  of  (31)  and  (32),  we  shall  obtain  a set  of  d-| 

and  a-|  by  substituting  the  displacement  components  of  a nodal  point 

into  the  left  hand  sides  of  (31)  and  (32).  Numerical  results  thus 

obtained  for  K-j  and  K2  vary  considerably  depending  on  the  locations  of 

nodal  points  and  the  values  of  displacement  components  used.  From 

(28)  and  (29),  the  displacement  components  of  a collapsed  triangular 

1/2  3/2 

element  are  functions  of  r , r and  r . If  u*  and  v*  designate  the 
leading  portion  of  displacements  (r^  term  only)  in  (28)  and  (29), 
we  have  specifically 
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(36) 


i 1/2  . v 

u*(n  = -1)  = g (r/S,)  (18u2  - 9u3  + 2u4)  | 

u*(n  = 1)  = 1 ( r/A ) 1 ' 2(18ug  - 9u8  + 2u7)  j 

and  similar  expressions  for  v*  (replacing  u by  v).  The  stress 
intensity  factors  obtained  from  (35)  with  d-|  and  a-j  computed  from 
(31)  and  (32)  using  u*,  v*  on  the  left  hand  sides  and  using  only  the 
v^r  term  in  the  right  hand  sides,  are  independent  of  r.  The  use  of 
the  leading  portion  of  displacements  u*  or  v*  on  the  left  hand  sides 
of  (31)  and  (32)  is  suggested  by  Tracey  [14]  and  discussed  by 
Barsoum  [15].  For  a mode  I (or  mode  II)  crack,  all  a's  (or  d's)  of 
(31)  and  (32)  vanish.  The  stress  intensity  factor  K-|  (or  K2)  can  be 
obtained  by  either  (31)  or  (32). 

2.  Two-Term  Expansion: 

1 12 

In  (31)  and  (32),  if  r and  r terms  are  considered,  we  have 

four  unknown  constants  d-| , a-j , d^,  a2  to  be  determined.  Four  displace- 
ment components  of  any  two  nodal  points  should  suffice  to  compute  K-, 
and  K2.  In  practice,  we  use  either  the  pair  of  nodes  (2,9)  or  (3,8) 
of  Figure  2.  Two  nodes  of  different  r such  as  (2,3)  or  (8,9)  usually 
give  poorer  results.  On  the  left  hand  side  of  (28)  and  (29),  we  may 
use  actual  nodal  displacements  (u2>v2),  (u9,vg)  or  we  may  use  (u2**, 
v2**)  and  (ug**,Vg**)  where  u2**  is  the  part  of  u2  corresponding  to 
r and  r terms  in  (28).  For  n = ± 1 

^Tracey,  D.  M.,  '’Discussion  of  'On  the  Use  of  Isoparametric  Finite 
Elements  in  Linear  Fracture  Mechanics'  by  R.  S.  Barsoum,"  Int. 

Journal  for  Numerical  Methods  in  Engineering,  Vol . 11,  1977. 

^Barsoum,  R.  S.,  "Author's  Reply  to  the  Discussion  by  Tracey,"  Int. 
Journal  for  Numerical  Methods  in  Engineering,  Vol.  11,  1977. 
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11/2  9 

u**(n  = -1)  = j (r/Z)  (18u2  - 9u3  + 2u4)  - j (r/fc)(5u2  - 4u3  + u4) 

11/2  9 

u**(n  = 1)  = 2 (r/Z)  (18ug  - 9u8  + 2uy)  - j (r/Z)(5ug  - 4uQ  + uy) 

(37) 

For  a mode  I crack,  both  a-|  and  a^  vanish.  Only  two  displacement 
components  of  a node  are  needed  to  determine  d-j  and  d„.  Any  of  the 
following  pairs  may  be  used  for  this  purpose:  (u2>v2),  (u3>v3),  (u^, 
Vg) , (u8,v8). 

3.  Four  Term  Expansion: 

If  we  take  r’^  up  to  r2  terms  in  (31)  and  (32),  the  eight 
constants  d^  and  a^ , i = 1,  2,  3,  4 are  to  be  determined  from  displace- 
ment components  of  four  nodal  points.  We  may  take  two  mid-nodes  ( 1/9 
and  4£/9  from  the  tip)  of  the  two  sides  of  a collapsed  triangular 
element  as  the  four  nodal  points,  or  we  may  take  four  consecutive 
nodes  of  the  same  radius  from  the  tip  (r  = £/9  or  r = 4£/9).  For  the 
three  collapsed  triangular  elements  around  a mode  I crack  shown  in 
Figure  5(a),  displacement  components  are  taken  from  one  of  the 
following  three  groups:  (i)  nodes  11,  12,  13,  14;  (ii)  nodes  15,  16, 
17,  18;  (iii)  nodes  11,  12,  15,  16  of  the  element  (1),  nodes  12,  13, 

16,  17  of  the  element  (2)  and  nodes  13,  14,  17,  18  of  the  element  (3). 

4 . Collocation  Method: 

Let  u_,  vs  be  displacements  from  the  asymptotic  expansions  (31) 
and  (32)  and  let  u and  v be  displacements  from  finite  elements  given 
by  (28)  and  (29).  For  an  arbitrarily  fixed  r,  we  define 
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e = z {[us(01)  - u^)]2  + Cvs (ei ) - ve(0i)]2} 


(38) 


The  unknown  coefficients  d's  and  a's  in  (31)  and  (32)  are  the  set 


which  minimizes  e.  In  other  words,  d's  and  a's  are  solutions  of 

\ 


3e  = 
3d,* 


= 0 
3a  n* 


N 


> j = 1,2,...,  n < j 


(39) 


where  n = 2p,  p is  the  highest  exponent  of  r in  the  asymptotic 
expansion. 


NASTRAN  IMPLEMENTATION 

The  NASTRAN  implementation  of  the  12-node  quadrilateral  follows 
that  of  the  8-node  quadrilateral  as  described  in  [10].  The  dummy 
users  element  facility  of  NASTRAN  is  used.  This  requires  coding 
routines  to  calculate  element  stiffness  matrices  and  stress  recovery 
computations.  Modifications  to  existing  NASTRAN  source  codes  are 
made  to  provide  proper  output  formats  for  the  element.  Three-point 
Gaussian  quadrature  is  normally  used  to  calculate  each  partial 
integration  of  the  double  integral  (4).  As  an  option  a four-point 
Gaussian  quadrature  may  be  used  instead.  All  stiffness  computations 
are  performed  in  double  precision  while  stress  recovery  is  performed 
in  single  precision.  Element  stiffness  matrix  computation  requires 
10  seconds/element  on  an  IBM  360/44.  Stress  intensity  factors  are 
calculated  from  nodal  displacements  by  various  techniques.  Equations 


^Hussain,  M.  A.,  Lorensen,  W.  E.,  and  Pflegl,  G.,  "The  Quarter-Point 
Quadratic  Isoparametric  Element  As  a Singular  Element  for  Crack 
Problems,"  NASA  TM-X-3428,  1976,  p.  419. 
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(40)  are  finally  adopted  to  compute  stress  intensity  factors  for 
mode  I,  mode  II  and  mixed  mode  cracks. 

NUMERICAL  RESULTS 

Computer  program  APES  [16]  utilizes  the  same  12-node  quadrilateral 
isoparametric  elements  shown  in  Figure  1.  The  high  order  element 
greatly  reduces  the  total  number  of  elements  as  well  as  nodal  points 
needed  to  model  either  elasticity  or  fracture  problems.  The  program 
has  been  designed  primarily  for  user  convenience.  There  are  many 
convenient  features  such  as  the  automatic  generation  of  middle  nodes 
of  a side  which  is  a straight  line,  and  the  automatic  computation  of 
nodal  force  for  a given  distributed  and/or  concentrated  tractions. 

However,  APES  requires  the  use  of  special  crack  tip  elements  for 
fracture  problems.  There  are  two  different  types  of  crack  tip  elements 
used  in  APES.  One  is  a small  circular  core  element  which  completely 
encloses  a crack  tip  and  which  reproduces  the  singular  nature  of  the 
strains  there.  The  other  consists  of  several  "enriched"  12-node 
elements  around  a crack  tip  all  having  a corner  node  corresponding  to 
the  tip.  In  an  enriched  element,  the  displacement  assumption  is 
augmented  by  the  leading  terms  of  the  singular  near  field  solution 
(an  extension  of  the  work  of  Benzley  and  Beisinger  [7]).  Both  models 
lead  to  about  the  same  high  degree  of  accuracy  in  stress  intensity 
factors  K-|  and 

7Benzley,  S.  E.  and  Beisinger,  A.  E.,  "Chiles  - A Finite  Element 
Computer  Program  That  Calculates  the  Intensities  of  Linear  Elastic 
Singularities,"  Sandia  Laboratories,  Technical  Report  SLA-73-0894,  1973. 

^Gifford,  L.  N.,  "Apes  - Second  Generation  Two-Dimensional  Fracture 
Mechanics  and  Stress  Analysis  by  Finite  Elements,"  Naval  Ship  Research 
and  Development  Center,  Report  4799,  December  1975. 
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In  this  report,  we  have  used  the  collapsed  triangular  elements  to 
eliminate  the  use  of  these  special  crack  tip  elements.  The  same 
displacement  functions  are  used  for  quadrilateral  elements  and  trian- 
gular elements,  and  there  is  no  question  regarding  the  continuity  of 
displacement  between  the  special  crack  tip  elements  and  conventional 
QUAD-12  elements.  There  is  no  need  for  the  use  of  8 x 8 Gaussian 
quadrature  for  numerical  integration  of  element  stiffness  matrix  as  is 
required  for  the  enriched  QUAD-12  elements.  From  the  numerical  results, 
the  collapsed  triangular  elements  as  singular  crack  tip  elements  will 
be  assessed. 

Three  mode  I crack  problems  and  one  mixed  mode  crack  problem  are 
chosen  for  numerical  computation  of  stress  intensity  factors.  The 
geometries  and  loads  of  mode  I tension  test  specimens  are  given  in 
Figure  3.  The  idealization  of  a half  of  the  single  edge  crack  is 
shown  in  Figure  4.  Similar  idealization  is  used  also  for  a quadrant 
of  a center  crack  or  a double  edge  crack.  Three  collapsed  triangular 
elements  surrounding  a mode  I crack  tip  are  shown  in  Figure  5(a).  The 
same  idealization  is  used  for  NASTRAN  and  APES  (with  collapsed  trian- 
gular elements).  For  comparison,  a circular  core  with  10-core  nodes 
around  the  crack  tip  is  also  used  in  APES  and  is  shown  in  Figure  5(b). 
Around  a mixed  mode  crack  tip,  six  collapsed  triangular  elements, 
shown  in  Figure  6(a),  are  used  and  the  corresponding  singular  'core' 
element  surrounded  by  six  QUAD-12  elements  are  shown  in  Figure  6(b). 

A 45°  slant  edge  crack,  in  a rectangular  panel  under  tension,  is  taken 
as  an  example  of  mixed  mode  cracks  and  is  shown  in  Figure  7.  The 
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sixteen-element  idealization  of  the  cracked  panel  is  also  shown  in  the 
same  figure,  and  the  six  elements  around  the  crack  tip  are  enlarged 
and  shown  in  Figure  6. 

Our  first  goal  in  the  numerical  analysis  is  to  find  a simple  and 
accurate  way  to  estimate  the  stress  intensity  factor  from  nodal 
displacements  obtained  from  the  finite  element  method.  For  the 
purpose  of  comoarison,  stress  intensity  factors  obtained  from  the 
finite  element  method  are  normalized  by  corresponding  values  which  a^e 
considered  a exact.  Referring  to  Figure  3,  K-|  = 1.966  [17]  is  taken 
as  exact  for  the  central  crack  and  K-j  = 2.00  [18]  for  the  double  edge 
crack.  For  the  single  edge  crack  K-|  = 3.728  for  a/b  = 0.4,  H/b  = 4 
and  K-|  = 5.009  for  a/b  = 0.5,  H/b  = 3 (first  F(a/b)  on  page  2.11  of 
[19]).  Normalized  strec  . intensity  factors  for  the  three  mode  I crack 
specimens  are  computed  Ky  APES,  using  singular  'core'  el  erne:  t with 
10-core  nodes,  Figure  5(b).  Results  are  tabulated  in  Table  1 for 
h/rQ  = 6 and  rQ  = 0.01,  0.02  and  0,03.  The  overall  results  are  better 
for  rQ/a  = 0.01.  Here  and  in  tables  1 through  8,  r0  = 0.01,  0.015,... 
should  be  understood  as  r0/a  = 0.01,  0.015,...  since  a = 1 is  used  for 
all  mode  I crack  examples. 

l^Isida,  M.,  "Analysis  of  Stress  Intensity  Factors  for  the  Tension 
of  a Centrally  Cracked  Strip  with  Stiffened  Edges,"  Engineering 
Fracture  Mechanics,  Vol . 5,  1973. 

18Brown,  W.  F.,  and  Srawley,  J.  E.,  "Plane  Strain  Crack  Toughness 
Testing  of  High  Strength  Metallic  Materials,"  ASTM  STP-410,  1966. 

19Tada,  H.,  Paris,  P.,  and  Irwin,  G.,  The  Stress  Analysis  of  Cracks 
Handbook,  Del  Research  Corp.,  1973. 
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TABLE  1.  RATIOS  OF  K]  (APES)  TO  K]  (EXACT) 


r0  = 0.01 

r0  = 0.02 

r0  = 0.03 

Center  crack  (exact  K-|  = 1.966) 

0.996 

0.976 

0.969 

Double  edge  crack  (exact  K-j  = 2.00) 

1.001 

0.993 

0.989 

Single  edge  crack  (exact  K-|  = 5.009) 

0.974 

0.964 

0.966 

Replacing  the  special  core  element  and  the  regular  12-node,  quad- 
rilateral elements  surrounding  the  core  with  collapsed  triangular 
elements,  we  can  use  the  APES  program  for  general  structures  (no  cracks) 
to  obtain  displacement  components  at  all  nodes.  Stress  intensity 
factors  are  then  obtained  from  displacements  at  nodes  close  to  the 
crack  tip  by  various  techniques  mentioned  previously.  If  only  one 
term  expansion  is  used,  the  ratios  of  K-|  (Finite  Element)  to  K-|  (Exact) 
for  the  three  tension  test  specimens  are  given  in  Tables  2-4.  The  use 
of  v*  (leading  term  of  v)  does  not  always  give  better  results  (e.g. 

Table  4).  Using  two  term  expansion,  the  same  ratios  are  given  in 
Tables  5-7.  The  results  from  two  term  expansion  show  little  improve- 
ment. Similar  results  using  four  term  expansion  are  tabulated  in 
Table  8.  The  values  in  the  last  column  of  Table  8 are  the  linear 
average  of  three  values  of  K-j  obtained  from  three  different  elements. 
Careful  study  of  results  tabulated  in  Tables  2 through  8 reveals  that 
one  term  expansion  using  v at  node  14  or  18,  Figure  5(a),  with  rQ  = 

1%  - 2%  of  the  crack  length  is  the  simplest  way  to  estimate  K-j  and 
K-j  thus  obtained  is  quite  accurate.  This  technique  is  adopted  in  our 
NASTRAN  to  compute  K-| . 
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BY  TWO-TERM  EXPANSION  USING  APES 
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The  collocation  method  requires  the  computation  of  displacements 
along  r=r0  at  points  between  nodes.  For  the  single  edge  tension 
specimen  with  a/b=0.5,  H/b=3.0,  the  numerical  results  obtained  from 
the  collocation  method  are  given  in  Table  9.  It  can  be  seen  that 
values  of  K-j  remain  nearly  a constant  no  matter  how  many  terms  or  how 
many  mid-points  are  taken.  Because  of  more  computations  involved  and 
since  as  good  or  even  better  results  can  be  obtained  by  other  simpler 
methods,  the  collocation  method  has  not  been  pursued  further. 

TABLE  9.  K-|  (FINITE  ELEMENT)  FOR  SINGLE  EDGE  CRACK  USING 

COLLOCATION  METHOD.  NODAL  DISPLACEMENTS  OBTAINED  FROM  APES 
WITH  3-COLLAPSED  TRIANGULAR  ELEMENTS.  COLLOCATION  POINTS 
ARE  EQUALLY  SPACED  ON  r=0.01 


No.  of  Mid=Points 
Between  Nodes 

Four-Term  Expansion 

Eight-Term  Expansion 

General 

Mode  I (a's=0) 

General 

Mode  l(a^=0) 

0 

4.90392 

4.85093 

2 

•3 

4.82446 

4.82966 

o 

5 

4.82815 

4.83417 

4.81972 

4.82762 

9 

4.82753 

4.83352 

11 

4.82743 

4.83337 

14 

4.82737 

4.83322 

19 

4.82734 

4.83309 

Our  second  goal  in  the  numerical  computation  is  to  assess  the 
accuracy  of  NASTRAN  in  linear  fracture,  using  12-node,  isoparametric 
elements  with  collapsed  triangular  elements  around  a crack  tip  and  to 
examine  the  effect  of  multiple  constraints.  Normalized  stress  intensity 
factors  thus  obtained  are  given  in  Table  10.  It  indicates  a high  degree 
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of  accuracy,  and  the  effect  of  multiple  constraints  is  insignificant. 
The  multiple  constraints  tend  to  increase  the  stress  intensity  factors 
slightly.  In  Table  10,  by  "No"  multiple  constraint  we  mean  nodes  2 to 
10,  Figure  5(a),  are  free  to  move  in  both  the  local  x-  and  y-directions 
with  respect  to  node  1,  and  node  1 has  no  displacement  in  the  local 
y-di recti  on  due  to  the  symmetry  of  the  problem  (mode  I crack).  The 
column  "Yes"  gives  results  obtained  by  assuming  v-j  = v2  = •••  ■ v-|Q  = 0 
and  u-|  = u2  = •••  = u1Q. 

TABLE  10.  K1  (NASTRAN)/«i  (EXACT) 


ro/a 

0.01 

0.015 

0.02 

Multiple  Constraint 

No 

Yes 

No 

Yes 

No 

Yes 

Center  Crack 
a/b=0.4,  H/b=4.0 
Exact  K-|  =1.966 

0.981 

1.013 

0.982 

0.999 

0.983 

0.994 

Double  Edge  Crack 
a/b=0.4,  H/b=4.0 
Exact  K-j  =2 . 00 

1.000 

1.021 

0.998 

1 .007 

0.999 

1 .002 

Single  Edge  Crack 
a/b=0.4,  H/b=4.0 
Exact  K-j  =3 . 728 

0.980 

1.003 

0.980 

0.991 

0.982 

0.988 

Another  goal  in  the  numerical  computation  is  to  compare  results 
using  the  concept  of  collapsed  triangular  elements  in  APES  with  results 
of  APES  with  singular  'core'  element  around  the  crack  tip.  Ratios  of 
K-|  (APES)  to  K,  (EXACT)  are  given  in  Table  11  for  the  tension  test 
specimens  with  rQ/a  = 0.01.  Results  in  the  first  column  of  Table  11 
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are  obtained  by  using  three  collapsed  triangular  Elements  around  the 
crack  tip  where  nodes  1 to  10,  Figure  5(a),  coincide.  All  vertical 
displacements  of  nodes  1 to  10  are  assumed  to  vanish  while  they  are 
free  to  slide  with  respect  to  one  another  in  the  x-direction.  Results 
in  the  second  column  are  computed  by  using  a singular  'core'  element 
surrounded  by  three  quadrilateral  elements  as  shown  in  Figure  5(b) 
with  h/rQ  = 6.  fable  11  shows  values  of  stress  intensity  factors, 
obtained  by  collapsed  triangular  element,  are  as  accurate  as  those 
obtained  by  using  singular  'core'  element. 


TABLE  11.  «i  (APES)/«i  (EXACT) 


Collapsed  Triangular 
Elements,  Fiq.  5(a) 

Singular  "Core" 
Elements,  Fiq.  5(b) 

Center  Crack 
a/b=0.4,  H/b=4.0 
Exact  l<i  =1.966 

0.998 

0.996 

Double  Edge  Crack 
a/b=0.4,  H/b=4.0 
Exact  K-|=2.00 

1.005 

1.011 

Single  Edge  Crack 
a/b=0.5,  H/b=3.0 
Exact  K-|=5.009 

0.972 

0.974 

In  NASTRAN,  it  is  the  user's  choice  to  apply  the  multiple  constraint 
conditions,  equations  (27),  at  the  crack  tip.  But  in  APES,  it  is  not 
yet  available  for  the  application  of  multiple  constraints.  Since  the 
effect  of  multiple  constraints  is  so  small,  Table  10,  that  results 
from  APES  with  collapsed  triangular  elements  are  considered  correct 
even  if  the  multiple  constraint  conditions  (27)  are  not  satisfied. 
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For  the  integration  of  the  element  stiffness  matrices,  a 3 x 3 
Gaussian  quadrature  is  used  in  all  NASTRAN  results  in  this  report  while 
a 4 x 4 Gaussian  quadrature  is  used  in  APES. 

For  a mixed  mode  crack,  six  collapsed  triangular  elements  shown 
in  Figure  6(a)  are  used.  An  effective  way  to  estimate  the  stress 
intensity  factors  is  to  use  v(r0,-Tr)  for  K-j  and  u(r0,Tr)  and  u(r0,-ir) 
for  K2. 


K-|  (tt) 
K2(tt) 


/2tt  2G  v(ro»ir) 
(k+1) 

/2tt  2G  u(r0,ir) 


Ki  = j (K-|  (tt)  + K-|  (-tt)) 


-/2tt  2G  v(r0,-7r) 

K-i(-tt)  = — 7—77 

1 /Tq  (K+1)  \ 


K2(-7t) 


K2  = 2 


-y^TT  2G  u(r0,-7T) 
v^O  (*+1) 
(K2(tt)  + K2(-tt)) 


) (40) 


J 


where  u(r0,7r)  and  v(rQ,Tr)  are  displacement  components  of  node  26 
relative  to  node  19  in  the  direction  of  parallel  and  normal  to  the 
crack  face;  rQ  is  the  distance  between  nodes  26  and  19.  u(r0,-7r)  and 

v(r0,-7r)  are  the  same  but  of  node  20  relative  to  node  1. 

For  a 45°  edge  crack  shown  in  Figure  7,  NASTRAN  results  for  K-| 
and  K2  are  tabulated  in  Table  12  for  rQ/a=0.01  and  for  various  other 
conditions.  Again  the  multiple  constraint  conditions,  namely  u-|  = U2  = 

•••  = u 1 g and  v-|  = v2  = •••  = v-jg,  give  little  effect  on  values  of  K-j 
and  «2.  An  obliqued  edge  crack  in  a rectangular  panel  under  uniform 
tension  is  solved  by  Freese  using  a modified  mapping  collocation  method 
[20].  K.|  and  K^,  read  from  Bowie's  graphs  (Figure  1.16a  and  1.16b  of 

20Bowie,  0.  L.,  "Solutions  of  Plane  Crack  Problems  by  Mapping  Technique," 
in  Mechanics  of  Fracture,  1,  Edited  by  G.  C.  Sih,  Noordhoff  International 
Publishing,  Leyden,  The  Netherlands,  1973. 
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[20]),  are  approximately  1.86  and  0.88  respectively.  Numerical  results 
of  the  same  problem  computed  by  APES,  using  special  crack  tip  elements 
and  various  idealizations,  are  given  in  Figure  12  of  [16]. 


TABLE  12.  K1  AND  K2  FOR  45°  EDGE  CRACK  BY  NASTRAN 


B.C. 

Integration 

Multiple  Constraint 

K1 

k2 

1 

3x3 

No 

1.89 

0.95 

1 

3x3 

Yes 

1.89 

0.96 

1 

4x4 

No 

1.83 

0.92 

2 

4x4 

Yes 

1.84 

0.93 

THE  STABILITY  OF  COLLAPSED  TRIANGULAR  ELEMENTS 


In  a recent  report  by  Hussain  and  Lorensen  [21],  it  was  found  that 
a slight  perturbation  in  placing  the  mid-side  node  opposite  to  the 
crack  tip  for  a collapsed  8-node,  quadrilateral  element  led  to  unstable 
results  in  stress  intensity  factors.  This  unstability  can  be  shown  in 
the  collapsed  12-node,  quadrilateral  element  if  one  or  both  middle 
nodes  of  the  side  opposite  to  the  crack  tip  are  slightly  perturbed  from 
their  nominal  positions. 

Let  node  5 be  perturbed  as  shown  in  Figure  8.  Denoting  the 
perturbed  quantities  with  an  asterix,  we  have 


^Gifford,  L.  N.,  "Apes  - Second  Generation  Two-Dimensional  Fracture 
Mechanics  and  Stress  Analysis  by  Finite  Elements,"  Naval  Ship  Research 
and  Development  Center,  Report  4799,  December  1975. 

20 

Bowie,  0.  L.,  "Solutions  of  Plane  Crack  Problems  by  Mapping  Technique,1 
in  Mechanics  of  Fracture,  1,  Edited  by  G.  C.  Sih,  Noordhoff  Inter- 
national Publishing,  Leyden,  The  Netherlands,  1973. 

21 

Hussain,  M.  A.  and  Lorensen,  W.  E.,  "Isoparametric  Elements  As 
Singular  Elements  for  Crack  Problems,"  Watervliet  Arsenal  Technical 
Report,  to  be  published. 
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(41) 


* 2 + cosa  . 
V*  = — r + E 


yj/t  - s' 

A general  point  (x,y)  given  by  equation  (18)  will  be  displaced  at 


x*/J l = ^ (l+5)2[(l-n)  + (l+n)cosa]  + e (l+5)(l-n2)(l-3n)  (42) 

y*/SL  = jjj-  (l+?)2(l+n)sina  + e'  ■Jr  (l+£)(l-n2)(l-3n)  (43) 

Along  the  line  n = -1/3,  and  replacing  y*  with  r sine  in  (40),  we  have 


1 + 5 


sina 


1 n r 4sin9sinou 

U l 3c12  1 


1/2 


(44) 


Since  (1H)  is  a common  factor  in  displacement  components,  equations 
(28)  and  (29),  it  is  seen  that  the  singularity  required,  for  the 
crack  problems  disappears  along  at  least  the  ray  n = -1/3  in  the 
collapsed  triangular  case. 

As  a numerical  example,  the  central  crack  tension  specimen  of 
Figure  3(a)  is  again  used.  If  the  idealization  remained  the  same  as 
shown  in  Figure  4,  except  that  the  collapsed  elements  of  Figure  5 were 
replaced  by  those  of  Figure  8(b)  (where  nodal  points  20,  21,  23,  24, 

26  and  27  are  on  a circular  arc),  the  computed  stress  intensity  factor 
changed  from  its  almost  exact  value  K-j  = 1.962  to  K-|  = 1.421  (nearly 
30%  error).  If  only  nodal  points  26  and  27  were  perturbed  to  their 
new  locations  in  Figure  8(b),  the  stress  intensity  factor  would  become 
K-|  = 1 .457  (a  26%  error) . 
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CONCLUSIONS 


By  a simple  manner,  the  12-node  isoparametric  elements  can  be 
used  to  form  a singular  element  for  two-dimensional,  elastic,  fracture 
mechanics  analysis.  The  elements  have  been  successfully  implemented 
in  NASTRAN,  which  can  now  be  more  efficiently  used  for  more  accurate 
prediction  of  stress  intensity  factors  of  complicated  crack  problems. 
The  middle  nodes  of  the  side  opposite  to  a crack  tip  in  a collapsed 
triangular  element  should  be  accurately  located  to  avoid  unstable 
results.  The  extension  of  collapsed  triangular  elements  as  singular 
elements  to  three-dimensional  brick  elements  can  be  easily  done  as  in 
[9,10]. 


Q 

Barsoum,  R.  S.,  "On  the  Use  of  Isoparametric  Finite  Elements  in  Linear 
Fracture  Mechanics,"  International  Journal  for  Numerical  Methods  in 
Engineering,  Vol . 10,  1976. 

10Hussain,  M.  A.,  Lorensen,  W.  E.,  and  Pflegl,  G.,  "The  Quarter-Point 
Quadratic  Isoparametric  Element  As  a Singular  Element  for  Crack 
Problems," NASA  TM-X-3428,  1976,  p.  419. 
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X 


N1  = J2  d-nD  (1-?)  [-10+9CC2+ri2)] 
N2  = J2  U-n)  (i-52)  (i-3?) 

N3  = ~k  d-nD  (1-^2)  (1+3^D 
N4  = J2  (i-n)  (i+5) [-10+9 C?2+n2) ] 
n5  = (i+5) (i-n2) (i-3n) 

n6  = J2  d+5)  (i-ti2)  (l+3n) 


Ns  = ^2  (i+n)d-52)(i+35) 

N9  = J2  d+n)  (i-52)  (i-35) 

Nio  = h Ci+n) Ci-?) [-io+9(52+n2j] 
Nn  = J2  Ci-?)  (i-n2)  d+3n) 
n12  = j2  (i-5) (i-n2) (i-3n) 


Figure  1.  Shape  Functions  and  Numbering  Sequence  For  a 12-Node 
Quadrilateral  Element. 
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Figure  2.  A Normalized  Square  in  (£>0)  Plane  Mapped  Into  a Collapsed 
Triangular  Element  in  (x,y)  Plane  with  the  side  £ = -1 
Degenerated  into  a Point  at  the  Crack  Tip. 


40 


t 

ro|<\j 


L 

^ICVj 


/ 0 


T 


(a)  Center  Crack 


I 0 


II 

'o 


-r  f 

f 

L_  . 

L H J 

H J 

(c)  Single-Edge  Crack 


Figure  3.  Three  Tension  Test  Specimens. 
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Figure  4 


Idealization  of  a Half  of  the  Single-Edge  Cracked 
Tension  Specimen. 
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Figure  5(a),  Three  Collapsed  Triangular  Elements  Surrounding 
a Mode  I Crack  Tip. 

(b) , Special  Core  Element  and  Three  Quadrilateral 
Elements  Surrounding  a Mode  I Crack  Tip. 
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Figure  6(a).  Six  Collapsed  Triangular  Elements  Surrounding  a 
Mixed  Mode  Crack  Tip. 

(b) . Special  Core  Element  and  Six  Quadrilateral  Elements 
Surrounding  a Mixed  Mode  Crack  Tip. 
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Figure  7.  Idealization  of  a 45-Degree  Slant  Edge  Cracked  Panel  in  Tension. 
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Figure  8 (a) . 
(b). 


Node  5 Perturbed  to  5*. 

Nodes  20,  21,  23,  24,  26,  27  Perturbed  From 
Their  Nominal  Positions. 
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US  ARMY  MATERIALS  6 MECHANICS 
RESEARCH  CENTER 

ATTN:,  TECH  LIB  - DRXMR-PL  2 

WATERTOWN,  MASS  02172 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARRADCOM,  ATTN:  BENET  WEAPONS  LABORATORY, 

DRDAR-LCB-TL,  WATERVLIET  ARSENAL,  WATERVLIET,  N.Y.  12189,  OF  ANY 
REQUIRED  CHANGES. 


